R Notebook

Authors:Martin Ostrowski, Deepa Varkey Version: 20180515

R4Oceanography

An introduction to R for plotting data on maps.

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

plot(cars)

cars is a pre-loaded dataset for demonstration purposes. To see what the dataset looks like execute ‘cars’ in a chunk

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

cars

   speed dist
1      4    2
2      4   10
3      7    4
4      7   22
5      8   16
6      9   10
7     10   18
8     10   26
9     10   34
10    11   17
11    11   28
12    12   14
13    12   20
14    12   24
15    12   28
16    13   26
17    13   34
18    13   34
19    13   46
20    14   26
21    14   36
22    14   60
23    14   80
24    15   20
25    15   26
26    15   54
27    16   32
28    16   40
29    17   32
30    17   40
31    17   50
32    18   42
33    18   56
34    18   76
35    18   84
36    19   36
37    19   46
38    19   68
39    20   32
40    20   48
41    20   52
42    20   56
43    20   64
44    22   66
45    23   54
46    24   70
47    24   92
48    24   93
49    24  120
50    25   85

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

---

Introduction to R

The aims of this practical are to develop core data handling, plotting and mapping skills using a programming lanaguage. R is the language that we will be using today. Other languages, such as ArcGIS, Matlab and Python, are also very capable of achieving the same, or better, results. Oh, by the way, R is free. Ultimately the choice of language will be yours. These notes are written in R Markdown format which will allow you to execute the code inline and produce a web friendly (.html) format preview.

The objectives of todays practical are to:

setup

• install the latest version of R on your own laptop (or one provided to you)
• install the latest version of RStudio
• Download the data package from github and move it to your Desktop folder

(20-30 min)

getting started

• Use R to plot an annotated map of Australia
• navigating the Rstudio environment
• loadng and manipulating data
• understanding basic data types: variables, vectors, dataframes, matrices etc.
• basic plots

• maps (45 min)

• exercises: practice customising plots

extending

• Obtain Satellite Sea Surface Temperataure and Chlorophyll data in the form of netcdf files
• Use R to plot SST on a map of Australia
• Use R to plot Chlorophyll concentration on a map

(30 min)

---

1. A guided introduction to R

(10 min)

2. Live coding, plot an annotated map of Australia

To plot data on a map you first need a map to work from. We can find the information by loading two packages, maps, and mapdata. Information in this section was sourced from Jeff Bowman

install.packages('maps')
install.packages('mapdata')
library(maps)
library(mapdata)

Let’s try plotting a map

map()   # low resolution map of the world

let’s limit the map region to Australia by addin limits for longitude x and latitude y using xlim=c(110,157), ylim=c(-48, -10)

map('worldHires', xlim=c(110,157), ylim=c(-48, -10), fill=T, col="grey")

map('worldHires', xlim=c(110,180), ylim=c(-48, -10), fill=T, col="grey", asp=1)

Now we can add data points.

Ocean Data

Australia operates an Integrated Marine Observing System which served data from ocean observations through the Australian Ocean Data Network

load data from the Integrated Marine Observing System National Reference Stations. IMOS operates 7 National Reference Stations around the coast of Australia.The coordinates for each station are given in the data file NRSlatlon.csv

Before loading the data we need to make sure that we are in the right directory with the setwd() command

setwd("~/Desktop/R4Oceanography/")

nrs.sites<-read.table("NRSlatlong.csv", h=T, sep=",")

check the data

nrs.sites

                      NRS code      lat     long
1            Port Hacking  PHB -34.1192 151.2267
2            Maria Island  MAI -42.5967 148.2333
3 North Stradbroke Island  NSI -27.3450 153.5620
4         Rottnest Island  ROT -32.0000 115.4167
5                 Yongala  YON -19.3085 147.6184
6         Kangaroo Island  KAI -35.8322 136.4473
7                  Darwin  DAR -12.4000 130.7680

Now try adding the data to the map using points() to add a plotting character pch= to an x,y coordinate

map('worldHires', xlim=c(110,180), ylim=c(-48, -10), fill=T, col="grey", asp=1)
points(nrs.sites$long, nrs.sites$lat)

These points are hard to see. We can make them bigger using the cex= option, and use a different plotting character (pch=16)

map('worldHires', xlim=c(110,180), ylim=c(-48, -10), fill=T, col="grey", asp=1)
points(nrs.sites$long, nrs.sites$lat, pch=16, cex=3)

This looks great.

Q how can we tell which site is which? Use text() to add text to an x,y position on a map

map ('worldHires', xlim=c(110,180), ylim=c(-48, -10), fill=T, col="grey", asp=1)
points (nrs.sites$long, nrs.sites$lat, pch=16, cex=5)

text (nrs.sites$long, nrs.sites$lat, nrs.sites$code, col="white")

Next, add the voyage stations from the RV Investigator voyage 04 in 2016

setwd("~/Desktop/R4Oceanography/")

The working directory was changed to /Users/mostrowski/Desktop/R4Oceanography inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

in16.sites<-read.table("IN16latlong.csv", h=T, sep=",")

Check the data

in16.sites

           CTD                  type    lat   long    col
1  IN16v04.001             Test cast -34.33 151.76    red
2  IN16v04.002         Eddy Transect -34.17 152.57    red
3  IN16v04.003         Eddy Transect -34.27 152.33    red
4  IN16v04.004         Eddy Transect -34.32 152.06    red
5  IN16v04.005         Eddy Transect -34.40 151.80    red
6  IN16v04.008   Process Stn Coastal -33.39 152.03    red
7  IN16v04.009   Process Stn Coastal -33.39 152.02    red
8  IN16v04.010                       -33.01 152.46   blue
9  IN16v04.011          Port Hacking -34.13 151.24  black
10 IN16v04.012            Upwelling? -36.32 150.28    red
11 IN16v04.013 Process Stn-  Coastal -36.25 150.30    red
12 IN16v04.014 Process Stn-  Coastal -36.26 150.29    red
13 IN16v04.015   Transect - Southern -36.23 151.27    red
14 IN16v04.016   Transect - Southern -36.25 151.93    red
15 IN16v04.017   Transect - Southern -36.25 153.26    red
16 IN16v04.018   Transect - Southern -36.25 154.25    red
17 IN16v04.019   Transect - Southern -37.00 154.00    red
18 IN16v04.020   Not the right water -37.00 152.99    red
19 IN16v04.021  Process Stn - Tasman -36.99 153.89    red
20 IN16v04.022  Process Stn - Tasman -36.99 153.89    red
21 IN16v04.023    Deep Cast - Tasman -36.89 152.55 purple
22 IN16v04.024                       -36.49 150.29    red
23 IN16v04.025                       -34.84 151.34    red
24 IN16v04.026     Process Stn - EAC -30.56 153.67    red
25 IN16v04.027     Process Stn - EAC -30.64 153.62    red
26 IN16v04.028     Drift Study - EAC -30.81 153.44 orange
27 IN16v04.029     Drift Study - EAC -30.97 153.39 orange
28 IN16v04.030     Drift Study - EAC -31.51 153.33 orange
29 IN16v04.031     Drift Study - EAC -31.65 153.32 orange
30 IN16v04.031     Drift Study - EAC -31.64 153.32 orange
31 IN16v04.032     Drift Study - EAC -31.84 153.32 orange
32 IN16v04.033     Drift Study - EAC -31.98 153.30 orange
33 IN16v04.034     Process Stn - EAC -32.46 153.18    red
34 IN16v04.035     Process Stn - EAC -32.48 153.16    red
35 IN16v04.035     Process Stn - EAC -32.47 153.17    red
36 IN16v04.036     Drift Study - EAC -32.57 153.12 orange
37 IN16v04.037     Drift Study - EAC -32.59 153.07 orange
38 IN16v04.038     Drift Study - EAC -32.76 152.99 orange
39 IN16v04.039     Drift Study - EAC -33.00 152.91 orange
40 IN16v04.040     Drift Study - EAC -33.05 152.88 orange
41 IN16v04.041  Process Stn - Tasman -32.47 154.49    red
42 IN16v04.042  Process Stn - Tasman -32.48 154.48    red
43 IN16v04.043     Transect - Middle -32.47 154.10    red
44 IN16v04.044     Transect - Middle -32.46 153.70    red
45 IN16v04.045     Transect - Middle -32.46 153.40    red
46 IN16v04.046     Transect - Middle -32.29 152.88    red
47 IN16v04.047     Process Stn - EAC -29.15 154.48    red
48 IN16v04.048     Process Stn - EAC -29.16 154.48    red
49 IN16v04.049   Transect - Northern -28.00 153.78    red
50 IN16v04.050   Transect - Northern -28.02 154.03    red
51 IN16v04.051   Transect - Northern -28.01 154.35    red
52 IN16v04.052   Transect - Northern -28.00 154.78    red
53 IN16v04.053   Transect - Northern -28.00 155.03    red
54 IN16v04.054             Deep Cast -27.92 155.15    red
55 IN16v04.055     Process Stn - EAC -28.00 155.04    red

Q: How many stations are there?

nrow(in16.sites)

[1] 55

Plot the stations on a map. Hint: change the longitude range.

map ('worldHires', xlim=c(145,160), ylim=c(-42, -20), fill=T, col="grey", asp=1)
title(xlab = 'Longitude',
      ylab = 'Latitude')

box()
grid()

points(in16.sites$long, in16.sites$lat, pch=16, col="red")

Exercises

1. Try changing the colour of the points using col= (“red”, “blue”, “black”)
2. Try changing the shape of the plotting characters using pch= (1-25)
3. Add the location of the nearset IMOS NRS
4. Colours are already defined in the column col (in16.sites$col). Q Can we use this to color the points?

map ('worldHires', xlim=c(145,160), ylim=c(-42, -20), fill=T, col="grey", asp=1)
points(in16.sites$long, in16.sites$lat, pch=16, col=in16.sites$col, cex=2)

3. Adding satellite data layers

Ocean Data Products specific to Australia are hosted on the AODN website. However, for this practical we will browse data from the NASA Ocean Color website, via the level 3 data browser. From the pulldown menus select Terra MODIS 4µ nightime, 4km resolution, 8 day composite. When hovering your mouse over the desired date you will see 3 download options in the bottom left, middle and right corners. Click on te leftee corner to download the .nc file and move it to your R4Oceanography directory on the desktop.

The .nc file extension signifies it is a netcdf file Network Common Data Form. Netcdf4 format is commonly used for high density arrays of data, such as those used in oceanography, climatology, taxonomy and GIS applications.

To enable R to work with netcdf files additional libraries need to be added

install.packages('ncdf4')
library(ncdf4)
install.packages('oce')
library(oce) #A Package for Oceanographic Analysis

What does a netcdf file look like?

sst01<-nc_open("~/Desktop/R4Oceanography/T20162572016264.L3m_8D_SST4_sst4_4km.nc")

sst01

File ~/Desktop/R4Oceanography/T20162572016264.L3m_8D_SST4_sst4_4km.nc (NC_FORMAT_NETCDF4):

     3 variables (excluding dimension variables):
        short sst4[lon,lat]   (Chunking: [1729,40])  (Compression: level 4)
            long_name: 4um Sea Surface Temperature
            scale_factor: 0.00499999988824129
            add_offset: 0
            units: degree_C
            standard_name: sea_surface_temperature
            _FillValue: -32767
            valid_min: -1000
            valid_max: 10000
            display_scale: linear
            display_min: -2
            display_max: 45
        byte qual_sst4[lon,lat]   (Chunking: [1729,40])  (Compression: level 4)
            long_name: Quality Levels, Sea Surface Temperature
            _FillValue: -1
            valid_min: 0
            valid_max: 5
        unsigned byte palette[eightbitcolor,rgb]   (Contiguous storage)  

     4 dimensions:
        lat  Size:4320
            long_name: Latitude
            units: degree_north
            _FillValue: -999
            valid_min: -90
            valid_max: 90
        lon  Size:8640
            long_name: Longitude
            units: degree_east
            _FillValue: -999
            valid_min: -180
            valid_max: 180
        rgb  Size:3
        eightbitcolor  Size:256

    65 global attributes:
        product_name: T20162572016264.L3m_8D_SST4_sst4_4km.nc
        instrument: MODIS
        title: HMODIST Level-3 Standard Mapped Image
        project: Ocean Biology Processing Group (NASA/GSFC/OBPG)
        platform: Terra
        temporal_range: 8-day
        processing_version: 2014.0
        date_created: 2016-10-06T10:24:39.000Z
        history: l3mapgen par=T20162572016264.L3m_8D_SST4_sst4_4km.nc.param 
        l2_flag_names: LAND,~HISOLZEN
        time_coverage_start: 2016-09-12T09:00:10.000Z
        time_coverage_end: 2016-09-20T11:45:08.000Z
        start_orbit_number: 89021
        end_orbit_number: 89139
        map_projection: Equidistant Cylindrical
        latitude_units: degrees_north
        longitude_units: degrees_east
        northernmost_latitude: 90
        southernmost_latitude: -90
        westernmost_longitude: -180
        easternmost_longitude: 180
        geospatial_lat_max: 90
        geospatial_lat_min: -90
        geospatial_lon_max: 180
        geospatial_lon_min: -180
        grid_mapping_name: latitude_longitude
        latitude_step: 0.0416666679084301
        longitude_step: 0.0416666679084301
        sw_point_latitude: -89.9791641235352
        sw_point_longitude: -179.97917175293
        geospatial_lon_resolution: 4.63831233978271
        geospatial_lat_resolution: 4.63831233978271
        geospatial_lat_units: degrees_north
        geospatial_lon_units: degrees_east
        spatialResolution: 4.64 km
        number_of_lines: 4320
        number_of_columns: 8640
        measure: Mean
        suggested_image_scaling_minimum: -2
        suggested_image_scaling_maximum: 45
        suggested_image_scaling_type: LINEAR
        suggested_image_scaling_applied: No
        _lastModified: 2016-10-06T10:24:39.000Z
        Conventions: CF-1.6
        institution: NASA Goddard Space Flight Center, Ocean Ecology Laboratory, Ocean Biology Processing Group
        standard_name_vocabulary: NetCDF Climate and Forecast (CF) Metadata Convention
        Metadata_Conventions: Unidata Dataset Discovery v1.0
        naming_authority: gov.nasa.gsfc.sci.oceandata
        id: T20162572016264.L3b_8D_SST4.nc/L3/T20162572016264.L3b_8D_SST4.nc
        license: http://science.nasa.gov/earth-science/earth-science-data/data-information-policy/
        creator_name: NASA/GSFC/OBPG
        publisher_name: NASA/GSFC/OBPG
        creator_email: data@oceancolor.gsfc.nasa.gov
        publisher_email: data@oceancolor.gsfc.nasa.gov
        creator_url: http://oceandata.sci.gsfc.nasa.gov
        publisher_url: http://oceandata.sci.gsfc.nasa.gov
        processing_level: L3 Mapped
        cdm_data_type: grid
        identifier_product_doi_authority: http://dx.doi.org
        identifier_product_doi: 10.5067/TERRA/MODIS_OC.2014.0
        keywords: Oceans > Ocean Temperature > Sea Surface Temperature
        keywords_vocabulary: NASA Global Change Master Directory (GCMD) Science Keywords
        data_bins: 15961783
        data_minimum: -1.43500006198883
        data_maximum: 35.131519317627

sst01.d<-ncvar_get(sst01, 'sst4');  ### check the variable name is correct, can be 'sst4', or other names
sst01.d[sst01.d>42.00072]<-NA  ### beware this value should be changed to NA in the matrix
sstlon <- ncvar_get(sst01, 'lon')
sstlat <- ncvar_get(sst01, 'lat')

An image with a colour pallette can be generated using imagep from the oce package

imagep(sstlon, sstlat, sst01.d, col=oceColorsTemperature(255), filledContour=T, missingColor=0, zlab='sst', xlim=c(140,165), ylim=c(-46, -25), zlim=c(10,28), asp=1);
map('worldHires', xlim=c(145,160), ylim=c(-46, -20), fill=T, col="lightgrey", add=T, border=NA)

title(xlab = 'Longitude',
      ylab = 'Latitude')

Next zoom in on the RV Investigator sites, and plot them on the map. BE patient.

imagep(sstlon, sstlat, sst01.d, col=oceColorsTemperature(255), filledContour=T, missingColor=0, zlab='sst (˚C)', xlim=c(145,160), ylim=c(-46, -25), zlim=c(10,28), asp=1)
points(in16.sites$long, in16.sites$lat, pch=16, col=in16.sites$col, cex=2)

map('worldHires', xlim=c(145,160), ylim=c(-46, -20), fill=T, col="lightgrey", add=T, border=NA)
title(xlab = 'Longitude',
      ylab = 'Latitude')

More exercises

1. find some satellite chlorophyll data, download it read into R and plot

That is it. Well Done!

Save the file. Preview the .html version. Keep this tutorial somewhere safe for future reference.

Where to get help and find useful material

• Attend a local Software Carpentry course run by Macquarie University
• Visit the Software Carpentry Website, find a course or access online lessons
• Google